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ABSTRACT 

A modified linear Kalman filter with deterministic 
forcing is used to improve the tracking capabilities of 
the Marine Corps AN/TPQ-27 remote tactical aircraft guidance 
system. Both sixth and ninth order filters are developed 
and used with the Precision and Coarse Guidance simulation 
programs. A technique for overcoming the effects of auto- 
pilot bias is presented and tested through simulation. 
Precision Guidance control is modified to utilize angle 
error and angle error rate to generate corrective commands. 
Significant improvement in aircraft state estimation, bomb- 
ing accuracy, and overall system response is shown. The 
Coarse Guidance algorithms and simulation program are nearly 
completely new and perform aircraft guidance to the bombing 
run with more than adequate precision under simulated condi- 
tions. The new version of the program is significantly less 
complex than the previous version and incorporates features 
which more realistically reflect actual conditions under 
which the system would be used. 
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I . 



INTRODUCTION 



A. AN/TPQ- 2 7 INTRODUCTION 

The AN/TPQ-27 system is a tactical aircraft guidance and 
control system used by the Marine Corps to guide strike force 
aircraft along a preplanned route, and then down a final leg 
for the purpose of performing precision bombing. 

The mission control and guidance is divided into two 
modes. The first mode, called Coarse Guidance, takes the 
aircraft from a TACAN entry point through a series of 
straight paths and command turns. These "legs" of the Coarse 
Guidance run are defined by the leg end points. The actual 
route to be followed is determined by tactical considerations. 
The last leg is also the bombing leg. After proper entrance 
onto the final approach to the bombing leg, the second mode 
of the mission guidance control takes command; this is called 
Precision Guidance. During the Precision Guidance mode, very 
accurate aircraft position estimates are developed to enable 
bomb placement calculations which will be accurate to errors 
on the order of feet, at the final impact point. The primary 
difference between these two modes, other than purpose, is 
the radar precision used to make the measurement of aircraft 
position. The resultant precision is also a function of the 
data rates of the respective radars, which differ signifi- 
cantly in this case. 

The initial efforts on this system have suffered from 
some problems which are also apparent from observation of the 
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results of software simulations of actual missions. The 
purpose of this study was to look into various techniques 
currently employed in the simulation programs, and attempt 
to improve the response and overall accuracy of the system 
through the use of different and/or improved algorithms. 

B. AN/TPQ -27 PROBLEM AREA DEFINITION 

The specific problem areas which were investigated were 
defined primarily through contact with personnel closely as- 
sociated with the system's performance. Other areas were 
noted as in need of improvement during the familiarization 
and simulation trial phases of study. 

1 . Problem Areas in Coarse Guidance 

Coarse guidance has suffered from many separate but 
related problems. Pilots have complained that the controls 
sent to the aircraft in the final moments before turning to 
a new leg have been violent. In addition, there have been 
complaints of not knowing where the aircraft was at any time 
other than having just exited from a "command turn," i.e., 
a turn from one leg onto a new leg, vice a course correction. 
The simulation program was found to be suffering from unneces- 
sary complexity in some areas, and was apparently in need of 
refined estimation and control procedures. 

2 . Problem Areas in Precision Guidance 

The single most prevalent complaint with the Precision 
Guidance program was the length of time required for the es- 
timates to stabilize in order than an accurate determination 
on exactly when the bomb should be dropped could be made. 
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Again, this and other related problems with Precision Guid- 
ance seemed due to inadequate aircraft estimation algorithms, 
and a control scheme which was too simplistic in design. 
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II. AIRCRAFT POSITION AND VELOCITY ESTIMATION 



A. BACKGROUND AND INTRODUCTION TO THE KALMAN FILTER 

The technique previously used to provide noise filtering 
on aircraft position and velocity was a standard alpha-beta 
filter with parameters chosen to yield the "optimal" track- 
ing capability in accordance with the theory developed in 
[1]. However, this reference also states that in the case 
when adaptive tracking is required, the a parameter should 
be permitted to vary with observed high frequency power 
fluctuations in the error signal 



e = x - x 

n pn 



CD 



where x is the state estimation, and x is the state pro - 
n ’ pn 

diction prior to measurement. Provisions for this variation 
were not included in the tracking algorithm which was imple- 
mented. In addition, the alpha-beta filter which was imple- 
mented was not an unbiased estimator of the aircraft state 
vector. This is due to the fact that as controls were used 
to cause changes in the free inertial model of the motion 
assumed by the alpha-beta filter, no corresponding changes 
were added to the filter states to account for this determin- 
istically added control. This accounts for the exceptionally 
large and prolonged transient errors which resulted from 
large control bank commands. 

The Kalman filter yields a minimum variance estimate of 
the state vector when the statistics of the noise are as 
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described below. This filter includes the effects of deter- 
ministic control commands to the aircraft to yield an esti- 
mator which is very nearly unbiased. The greatest improvement 
in estimation is yielded during the initial filter transient 
behavior. This is particularly critical in this application, 
in order to overcome the long filter settling period required 
by the alpha-beta filter. 

B. KALMAN FILTER ASSUMPTIONS AND GENERAL RECURSION EQUATIONS 
Application of the Kalman filter assumes that the discrete 
system under consideration satisfies 

X(k+1) = cj> (k)X(k) + W(k) (2) 

Z(k) = H (k) X (k) + V (k) (3) 

where X is an n x 1 state vector, Z is an m x 1 output vector, 
W is a zero-mean n x 1 vector of state excitation white noise, 
uncorrelated with the zero-mean additive white noise vector V, 
(j> is the state transition matrix (n x n) , and H is the m x n 
observation or measurement matrix. The assumed noise sta- 
tistics are 

E[V(k) V(j) T ] = R(k) 6(k,j) (4) 

E{r[W(k) W(j ) T ] r T } = Q (k) 6(k,j) . (5) 

E [V (k) W(j) T ] = 0 for all k,j (6) 

where 

6(k,j) - { 0 (7) 

1 k=j . 

The actual Kalman filter recursion equations are summarized 
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below, where X(k/j) denotes the estimate of the state X(k) 
based upon the j measurement observations Z (1) , Z (2) , . . . , Z ( j ) . 

P (k/k-1) = <|>(k,k-l)P(k-l/k-l)<Kk,k-l) T + Q(k) (8) 

G(k) = P (k/k - 1)H (k) [H (k) P (k/k-l)H (k) T + R(k)]' 1 (9) 
P (k/k) = P (k/k-1) - G(k)H(k)P(k/k-l) (10) 

X (k/k) = X(k/k-l) + G (k) [Z (k) - H (k) X (k/k- 1) ] (11) 

X(k/k-l) = <Kk,k-l)X(k-l/k-l) + T(k)U(k-l) (12) 

where P (k/k-1) is the covariance of error for the state pre- 

A 

diction vector X(k/k-l), P(k/k) is the covariance of error 

A 

matrix for the state estimation vector X(k/k), and G(k) is 
the gain matrix to be applied at the time of the k L measure- 
ment. Further detail on the summary and development of the 
Kalman filter equations is available in [2], [3], and [4], 

C. SELECTION OF REQUIRED SYSTEM MODEL ORDER 

The order of the filter refers to its capability to track 
a target exhibiting a particular type of motion without error, 
in a noiseless environment. For a single dimensional problem, 
a first order filter would estimate position only. Second 
order filters are capable of tracking a constant velocity 
target and estimating both position and velocity. Similarly, 
third order filters are capable of tracking a target exhibit- 
ing a constant acceleration profile, estimating position, 
velocity and acceleration in the process. 

Documentation provided on the original AN/TPQ-27 programs 
indicated that the aircraft would be flying a constant airspeed 
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profile.- This then indicates that the aircraft model might 
be appropriately chosen to be a free inertial (1/s 2 ) plant. 
Command orders to the aircraft in the form of bank angles 
serve to change the aircraft's heading in a totally deter- 
ministic manner, provided the true transfer function of the 
aircraft is known with respect to roll response. In view of 
the above, it was originally thought that it would be suffi- 
cient to implement second order Kalman filters to estimate 
position and velocity in each of the three coordinates, 
yielding a sixth order filter. 

It was later discovered that a problem known as auto- 
pilot bank angles bias exists with sufficient frequency and 
resultant imprecision that the order of the original feed- 
back loop was increased to compensate for the error through 
the use of a discrete integrator [5]. Concepts such as in- 
tegrators and digitally implemented lead- lag networks were 
to be avoided in the improved version of the simulation 
routines due to the resultant phase lags which they introduce. 
This is not to say that the steady-state response is in error, 
but simply that the time to reach that state is unsatisfac- 
torily long. When simulations were run with non- zero auto- 
pilot bias bank angles in conjunction with second order Kalman 
filters, significant errors resulted. 

A bias angle causes the aircraft to turn in a given direc- 
tion at a constant heading rate. It was postulated that it 
might be possible to estimate the bias angle, and induce an 
anti-bias in the estimator, but the noise on the bias estimate 
proved to be excessive. 
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Third order filters will estimate a constant acceleration. 
The second attempt to overcome the bias problem was to postu- 
late that the acceleration in the horizontal components would 
not change significantly over the relatively short periods of 
flight time in question. (Of course, if a bias existed and 
the aircraft was flown for a long enough time, the path flown 
would appear as a circle, and the third order filters could 
not possibly be satisfactory for that situation.) Thus, de- 
velopment of two separate filtering schemes was pursued. A 
sixth order (second order for each of three dimensions of 
motion) and a ninth order (third order for each of the three 
dimensions) filter were developed and tested. Use of the 
programs is very similar for each of the filters and is 
described below. 

D. DERIVATION OF THE COVARIANCE OF MEASUREMENT NOISE MATRIX, 

R(k) 

The Kalman Filter assumptions include linear relationships 
among measurements and states, as well as linear state tran- 
sition dynamics. The first of these is of concern at this 
point. The states of concern in measurement have been se- 
lected as Cartesian coordinates (x,y,z). However, the radar 
measures range, azimuth, and elevation. The assumed relation- 
ship between the states and the measured values are as shown 
in Fig. 1 and given by the equations below. 



x = R 


cos 0 


sin<j> 


(13) 


y = R 


cosG 


cosrj) 


(14) 


z = R 


s in0 




(15) 
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z 




x = R co s 0 ,sind> 

0 0 o T o 



* R o cos0 o cos * 0 
z o = R o sine o 



Figure 1. Illustration of assumed coordinate system used in 
the radar filters. 
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where R -is the range to the aircraft, 0 is the elevation 
angle, and <{> is the azimuth angle of the aircraft from 
North . 

Note that this form departs from that as given in (3) 
since the relationship between the measured variables and 
the states is a nonlinear one. If the states were directly 
observable, then (3) would appear as 



1 

rH 

N 




vl(k) 


z2(k) 


= H (k) X (k) + 


v2 (k) 


1 

N 

Od 

X* 

1 _ 




v3(k) 



where for a sixth order system 



(16) 



and 



H (k) 



1 

0 

0 



0 

0 

0 



0 

1 

0 



0 

0 

0 



0 

0 

1 



0 

0 

0 



(17) 



X(k) = [x(k) x(k) y (k) y(k) z (k) z (k) ] 



(18) 



Similarly, for a ninth order system 



X(k) = 



100000000 

000100000 

000000100 



and 



(19) 



( 20 ) 



X (k) = [x(k) x(k) x(k) y (k) y(k) y(k) z(k) z(k) 



Thus, in both cases the equations reduce to 



zl(k) 




x(k) 




"vl(k)‘ 


z 2 (k) 


= 


y(k) 




v2 (k) 


z 3 ( k ) 


' 


z(k) 

L J 




v 3 ( k ) 

— 



17 



z(k)] T . 



( 21 ) 



t h 

where vl"(k) is the k component of noise added to the x co- 

"til 

ordinate, v2 (k) is the x component of noise added to the 

"th. 

y coordinate, and v3(k) is the k component of noise added 
to the z coordinate. Due to the linearity of this problem, 
the three coordinate components might be considered as being 
statistically independent, in which case the R(k) matrix 
would probably be a diagonal array consisting of the indi- 
vidual coordinate measurement variances which would be con- 
stant for all k. For the nonlinear radar problem, the 
relationships for each k are 

zl = (R + n r ) cos(0 + n Q ) sinO + n^) (22) 

= x + vl 



(R + n r ) 


cos (6 + n e ) 005(4- + n^) 


t — \ 

to 


y + v2 






(R + n r ) 


sin (6 + n Q ) 


(24) 


z + v3 







where n , n n , and n, are incremental noise disturbances to the 
r ’ 0 ’ 4> 

true range, elevation and azimuth, respectively. Expanding 
(22) yields 

zl = (R + n ) [cos6 cosn 0 - sin6 sinn 0 ] • (25) 

[sin<+) cosn^ + cos4> sinn^] . 

It is assumed that n 0 and n^ are small angle perturbations, 
and therefore that 



cos n f] = cos n 
sin iig = n Q 



1 
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(26) 

(27a) 



sm n, = n, 

4> <f) 



(27b) 



(25) can then be written as 



zl = (R + n ) [cos0 sin<j) - n 0 sin6 sin4> 



(28) 



+ n^ cos0 cos«Ji - n^ n^ sin0 cost] 



= R cos0 sincf) + vl 



(29) 



where 



vl = 'RnQ sin0 sin4> + Rn^ cos0 cos<() - Rn^ n^ sin0 cos4> 

+ n r cos0 s-in<{> - n^ n^ sin0 sin4> + n r n^ cos0 cos<}> 



- n n n n, sin0 cos4>, 
r 0 r 



(30) 



Similar developments for the zl and z2 measurements yields 

v2 = -Rn n sin0 cos4> - Rn, cos0 sin* + Rn„ n, sin0 sin4> 

0 4> 0 <p 

+ n cos0 cos4> - n n„ sinu cos® - n n, cosO sin* 
r Y r 6 r 4> T 



+ n n Q n, sin0 sin4> 
r 0 4> 



(31) 



and 



v3 = Rn n cos0 + n sin0 + n n n sin0. 
0 r r 0 



(32) 



Repeating (4) , the equation for the covariance of measurement 
error matrix is 



R(k) = E < 



r 






vl(k) 


< 


v2 (k) 




v3(k)_ 







[ vl (k) v2(k) v3 (k) ] 



(33) 



J 

It is assumed that n^ is much smaller than R. Evaluating the 
diagonal terms of R(k) yields 



R(l,l) = R ? c 2 sin 2 0 sin ? 4> + R z o| cos 2 9 cos 2 4* 

+ R 2 0 q 0 2 sin 2 0 cos 2 4* + o 7 cos 2 0 sin 2 4> 



(34) 



19 



( 35 ) 



R(2;2) = R 2 a 2 sin 2 0 cos 2 <J> + R 2 cr| cos 2 0 sin 2 <}> 

+ Ra 2 a 2 sin 2 0 sin 2 <() + a 2 cos 2 0 cos 2 <f) 

0 <p r 

R(3 , 3) = R 2 a 2 cos 2 0 + a 2 sin 2 0. (36) 

/ 

The off-diagonal elements are simply the expected values of 
the cross product terms, and are computed in the same way 
yielding 



R(l,2) = R 2 a 2 (l - o 2 )(sin 2 0 sine}) cos4>) 

y y (37) 







+ (cr 2 


- R 2 a|)(cos 2 0 sin<{) cos4>) 




R( 2 


,3) 


= (a 2 - 
v r 


R 2 a 2 ) (sin0 cos0 


cos4>) 


(38) 


R(1 


,3) 


= (a 2 - 
v r 


R 2 Oq) ( sin0 cos0 


s in<h ) . 


(39) 


Due to 


the 


syiiime t ry 


of the R(k) arr 


’ ay , it is also 


-f~ -v-. -i -i -f- -t- 

L i C-4 ^ UiiU U 


R(2 




= R(l > 2) 






(40) 


R(3 


,d 


= R(l,3) 


and 




(41) 


R(3 


,2) 


= R(2,3) 






(42) 



Note that the R(k) matrix is in fact not a constant array, 
but is one which is state dependent, or rather R, 0, and <|> 
dependent . 

E. STATE PREDICTION EQUATIONS 

The state prediction equations are used to predict ahead 
from the current estimate to some arbitrary future point in 
time. Normally this time is that of the next measurement, 
however, this is not always the case. As a relevant example, 
in Coarse Guidance, it is required to predict ahead several 
times between radar samples, due to the long sampling interval 



2 0 



of the radar, and the need to precisely determine times to 
order command turns. 

1 . State Prediction in the Linear Case 

The analysis which follows will be addressed primarily 
to the ninth order case. The state vector is given in (20). 
The strictly linear prediction equations are given by (12) 
where 



<f>(k,k-l) = 



F (k) = 



1 


T 


ro|H 

ro 


1 0 


0 


0 


1 


0 


0 


0 




0 


1 


T 


1 0 


0 


0 


1 


0 


0 


0 


0 


0 


1 


1 0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


1 1 


T 


T 2 

2 


1 


0 


0 


0 


0 


0 


0 


I 0 


1 


T 


1 


0 


0 


0 


0 


0 


0 


1 0 


0 


i 


1 


Q 


Q 


0 


0 


0 


0 


r ' 
0 

1 


0 


0 


I 

1 


1 


T 


T 2 
2 


0 


0 


0 


1 

0 

1 


0 


0 


1 

1 


0 


1 


T 


0 


0 


0 


1 0 


0 


0 


1 

1 


0 


0 


1 






r 

T 3 


/ 6 




0 








0 








T 2 


/ 2 




0 








0 








T 


- 




0 








0 








0 






T 3 


/ 6 






0 








0 






T 2 


/ 2 






0 








0 






T 








0 








0 






0 








T 3 


/ 6 






0 






0 








T 2 


/ 2 






0 






0 








T 







(43) 



(44) 
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and 



U (k) = U x (k) a y (k) a z (k)] T . (45) 

Note that although in general the <}> and T matrices are func- 
tions of k, in this case they are not. This is an arbitrary 
choice. Normally, T represents the interval between sampling 
points and is a constant. If one wished to predict ahead by 
2T, he could either perform the prediction operation twice 
in succession, or simply compute a new cj> and T matrix using 
2T in place of T. The former technique is used in this study. 

The U(k) array represents the deterministic forcing 
for each of the three dimensions; the units are distance per 
sec 3 since the forcing is an acceleration rate. (For the 
sixth order filter U(k) is an acceleration (distance/sec 2 ).) 
This must be a function of k since at each time point the 
control to the aircraft may be different, and in general will 
be different. 

The above linear prediction equations of motion are 
those which would normally be used in a Kalman Filter. How- 
ever, it was found that these equations did not yield results 
of sufficient accuracy. A nonlinear technique was required 
to obtain' precise results. This is derived below, along with 
the equations for aircraft response to a bank angle command. 

2 . Aircraft Response to Bank Commands 

The aircraft is sent bank angle commands in order to 
cause heading changes. The flight profile is assumed to be 
a coordinated turn, in which there is no motion in the verti- 
cal plane, and the heading rate change is proportional to the 
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( 46 ) 



bank angle, as described in [6]. The relationship is 

¥ = (g/v) <f> a 

where \p is the heading rate in degrees/sec, g is the earth's 
gravitational constant (32.2 ft/sec 2 ) , V is the airspeed in 
ft/sec, and <f> is the actual angle' of bank. 

3 . 

The roll transfer function of the aircraft can be 
approximated by 

= V (ST b + D (47) 

/ 

as given in [6], where is the response time constant. 
is a function of the particular aircraft in use. The dis- 
cretized version of the solution to the differential equation 
resulting from (47) is 

- T / T b - T/T b 

*00 - * c (k-l)Cl - e D ) + *(k-l)e D (48) 

where T is the interval between predictions, <J> is the com- 
manded bank angle, and <p (k) is the actual bank angle at time 
T (k) . Then the turning rate is 

i(k) = (g/V) <J> (k) . (49) 



To get the incremental heading angle change over the time T, 
integrate (49) . 



AiHk) = 



* clt = (g/V) (p (k) dt 






= (g/V) 



i 
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T v 



(50) 
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$ c (k)t + (4>(k-l) - <j> c (k)) ^ r 
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Then 



= (g/V) [4> c (k)T 



( 50 ) 



(4>(k- 1) - 4» c (k))(T b )(l - e 




)]. 



i|> (k) = ^ (k- 1) + A<Kk) • (51) 

3 . State Prediction in the Nonlinear Case 

A ninth order filter is assumed for this analysis. 
The analysis includes effects of wind motion in x and y, and 
assumes knowledge of the wind components. Since the filter 
assumes a constant acceleration track, the acceleration pre- 
diction equations are 



x (k/k-1) 


= x(k-l/k-l) 


(52a) 


y (k/k-l) 


= y (k-l/k- 1) 


(52b) 


z (k/k-1) 


= z (k-l/k-1) . 


(52c) 



Since the new heading is known in terms of the command bank 
angle, the predicted velocities are 

A 

x (k/k- 1) = V (k- 1) s in (k) ] + W + i(k-l/k-l)T (53a) 
y (k/k-1) = V(k-l)cos[if»(k)] + w y + ?(k-l/k-l)T (53b) 
z (k/k-1) = z (k- 1/k- 1) T (53c) 



A 

where V(k-l) is the estimated air speed at time T(k-l), and 
W x and W are the estimated wind components. The predicted 
position is approximated through the use of numerical integra- 
tion using the Euler-Maclaurin summation formula, retaining 
only the first correction term. In general, this formula is 
given by [7] as 




f(t)dt : t J 0 
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The equations involving the integrals are 



x (k/k-1) 



y (k/k-1) 



z (k/k-1) 



x(k-l/k-l) + 



y (k-l/k-1) + 



z (k- 1/k- 1) + 



r 



x(t) dt 



T(k-l) 



p) _ 

J y(t) dt 

T(k-l) 

T(k) 

J z(t) dt 

T(k-l) 



(55a) 



(55b) 



(55c) 



Applying (54) yields 



x (k/k-1) 



= x (k-l/k-1) + | [x(k-l/k-l) + x (k/k-1) ] 

m 2 a a (56a) 

- jj [x t (k/k-1) - x t (k-l/k-1)] 



y (k/k-1) = y (k-l/k-1) + [y(k-l/k-l) + y(k/k-l)] 

m2 a a (56b) 

- {2 ty t (k/k-1) - y t (k-l/k-1)] 

A 

2(k/k-l) = T z (k/k-1) ' (56c) 

where x^. and y represent total accelerations in x and y. The 
simplification of the equation used to predict altitude is a 
result of the fact that there is no deterministic forcing in 
that direction. 

A discussion of the relationship between the accelera- 

A 

tion estimate X(k/k-l) and the total accelerations shown in 
(56a) and (56b) is required at this point. The purpose of 
estimating an acceleration is due to the fact that bias bank 
angles exist which tend to cause a continuous turning motion. 
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and thus, additional and unknown accelerations in the x and y 
directions. The prediction method as described above will 
place the aircraft at exactly the correct point with the 
correct velocity if there is no bias. In this case, the es- 
timator will estimate that zero additional acceleration is 
present. In the case in which an unknown bias bank angle is 
present, the prediction equations will not be moving the air- 
craft heading angle the correct number of degrees/sec and 
thus lag will develop in the filter. This lag will cause a 
non-zero acceleration to be estimated, and this value will 
subsequently be added appropriately to the respective vel- 
ocity components in order to compensate for the insufficient 
motion . 

In summary, the acceleration estimate is not an esti- 
mate of the total acceleration, but is only an estimate of 
that acceleration which is unknown, due possibly to a variety 
of causes, but due primarily to an unknown bias angle in the 
autopilot roll reference. Thus the derivatives with respect 
to time of the velocity components in (56a) and (56b) do not 
equate to the acceleration estimates. Rather, they may be 
computed by differentiating equations (53) with respect to 

A 

time and assuming that X(k-l/k-l) is approximately constant. 
This yields 

5 (k/k-1) = V(k-l) cos[y(k)]i(k) (57a) 

? t (k/k-l) = -V(k-l) sin[<Kk)]$i(k) (57b) 

and 
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(58a) 



'x (k-l/k-1) = V(k-l) cos (k- 1) ] #(k-l) 

(k-l/k-1) = -V(k-l) sin (k- 1) ] i(k-l). (58b) 

F. KALMAN FILTER IMPLEMENTATION 

Two separate Kalman filter subroutines were developed to 
simulate the operation and filtering of the radar processor. 
RADAR6 is the sixth order filter and RADAR9 is the ninth order 
filter. Their implementation in software is very similar; the 
differences are described below along with specific character- 
istics of the Kalman filter which apply to both versions. A 
general flow diagram for the filter is presented as Fig. 2. 

1 . Differences in RADAR6 and RADAR9 

The obvious difference between the two filters is 
that one is capable of tracking a maneuvering aircraft with 
nonzero autopilot bank angle bias, and the other is not. This 
is due to the fact that the sixth order filter does not include 
the acceleration states in x, y, and z. The price paid for 
this additional estimation capability is an increase in com- 
putation time and program size. As implemented, the ninth 
order filter requires about 32 percent more storage alloca- 
tion in memory than does RADAR6 ; RADAR6 requires about 76 K 
bytes of storage on an IBM 360-67. The increased computation- 
al time is difficult to judge since all timing data refers to 
the running of the overall program. The estimate for in- 
creased run time for the ninth order filter is on the order 
of 20 to 50 percent, depending on overall program run time. 

The increased run time and storage requirements result 
from the requirement that the ninth order filter must have 
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I igurc 2. General block diagram of Radar Subroutines. 
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arrays which are 9x9, while the sixth order filter requires 
only arrays which are 6x6; 18 of these arrays are involved 
resulting in a factor of 2.25 increase in array storage alone. 

The prediction calculations for the sixth order fil- 
ter differ from the ninth order filter only in the fact that 
the acceleration terms are missing in augmenting the pre- 
dicted velocity components. 

RADAR6 was used to develop the Coarse Guidance pro- 
gram since no mention of autopilot bias problems was found 
in the available documentation on that original system, [8] 
and [9]. RADAR9 was developed directly in response to the 
bias problem and was therefore used exclusively in the im- 
provement of the Precision Guidance simulation program. The 
two subroutines could be interchanged to work with the other 
simulation mode's program in a matter of minutes, should this 
be desired. 

2 . Initialization and Constant Array Calculations 

Upon entering either RADAR subroutine for the first 
time, logic passes control through a section of constant 
array calculation and definition of program constants. The 
constants specified at this point are primarily those which 
are used in bank and turn angle equations. These are a 
function of the aircraft type and update intervals and re- 
main constant throughout a given simulation run. 

A total of five arrays are defined at the program 
start and remain constant. An Identity array of order equal 
to the filter order is set up. The Measurement matrix, (II), 
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as given in either (17) or (19), depending on filter order is 
then defined. The state transition array, <f> , is considered 
constant for this application and is defined in accordance 
with the filter order. The Q matrix, which is a measure of 
the expected unknown random forcing to be applied to the sys- 
tem is computed using the T array as given in (44) (for the 
ninth order filter) in conjunction with the expected value 
of the random forcing array, W. For this study, W has been 
set to 0. The use and effects of non-zero values in this ar- 
ray are discussed below. 

The covariance of the initial state prediction vec- 
tor is set to 10 6 , a somewhat arbitrarily large number. The 
effect of choosing such a large variance on the states is to 
cause the filter to set the initial state estimation vector 
equal to the observation. In other words, the filter essen- 
tially ignores the a -priori information set in as initial 
conditions. This is a very typical method of initializing a 
linear Kalman filter, when little confidence is placed on 
any initial estimate of the states. The original simulation 
programs followed this practice (of zero initialization) and 
it was continued in the improved version. However, consider- 
able improvement in initial filter settling for the Kalman 
filters implemented could be achieved through the use of 
good initial conditions. The improved version of Coarse Guid- 
ance yields good estimates of the aircraft position, and in 
a "pass -off" to the Precision Guidance mode in a "live run" 
the states active in Coarse Guidance at the time would make 
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valid initial conditions for the Precision Guidance radar 
filter. 

G. PREDICTION ERROR, THE GAIN SCHEDULE, AND COVARIANCE OF 

AIRCRAFT MANEUVER 

The gain matrix, G(k), determines to what extent the 
data will be permitted to affect the state estimation. Dur- 
ing the initial few estimations, the data usually always 
plays a large role in determining states, due to the large 
uncertainties in the actual states specified by the large 
covariances of prediction and estimation. As more data is 
taken, the covariances fall, and there is less requirement 
to accept noisy data into the filter to update the state es- 
timate. Figure 3 shows a plot of the typical gain schedule 
for the x coordinate, C(l,l) term as a function of time. 
Values for the primary gain terms in y and z are similar. 

The curve is approximately exponential in shape, starts es- 
sentially at 1.0 at time zero, and decreases monotonically. 

As can be seen from (9) , the exact values in the gain matrix 
depend on a multitude of variables, including the states 
themselves in the form of the R matrix. 

Absolute compliance with the Kalman filter assumptions 
results in a gain schedule which is optimal in the sense that 
the state estimates \tfill be of minimum variance. The use of 
changing gains with time is implicit. Note that in the 
original simulation programs, only constant gain filters 
were used, or in a 'few instances, filters which started at 
one constant gain value and then switched to one different 
set of gains. 
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Figure 3. Illustration of an actual gain schedule decreasing 
as a function of time. 
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Time (seconds) 



The .gains are related to the filter bandwidth. High 
gains correspond to wide bandwidth since they "let in" nearly 
all of the measured data, including the noise on the data. 
Lower gains correspond to narrow band filters since very 
little data gets into the state estimation calculations. 

Use of low gains also resembles a narrow band low pass fil- 
ter in the phase lag which results in the state estimation 
vector due to an abrupt change in the actual states. When 
this results, large differences begin to develop in the pre- 
diction residual of the estimation equation (11) : 



E (k) = Z (k) - H(k)X(k/k-l) . (59) 

In this case, (59) simplifies to 

E x = x data • x t k/k -U (60a) 

E y * ^ata ' C60b 5 

E z z data - (60c) 



As the differences in equations (60) begin to be biased either 
positively or negatively over a period of time, the filter 
will begin to integrate to a new trajectory in state space to 
compensate for the fact. The lower the gains, the longer this 
process will take. 

If no random excitation noise, or no unknown forces or 
uncertainties are present in the system whose states are to 
be estimated, then the proper setting for W, the covariance 
of random state excitation, is 0. Since the Q array in (9) 
is 

Q (k) - r E [W W T jr r , (61) 
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if W is 0, then Q will also be 0. Examination of (9) in this 
case will reveal that as t->°° the gain schedule will go to zero, 
and as a result data will have little effect on the state es- 
timation process after a relatively short period of time after 
initialization. This is fine if all factors are known. 

In real systems, all factors affecting the states are 
rarely known. Uncertainties in AN/TPQ-27 might include wind 
velocity, exact aircraft roll response, biases in measurement 
equipment, and time variations of all of these. 

The goal in establishing improved simulation programs was 
ito devise techniques which would perform the required estima- 
tion and tracking assuming the above uncertainties did not 
exist to a significant extent. This attitude seemed to follow 
.that used in the generation of the original program. It was 
recognized that use of a gain schedule which goes to zero is 
not a likely useable solution to the problem; however, it did 
seem reasonable to try to generate a solution which worked with 
zero Q when the uncertainties did not exist. With this ac- 
complished, actual utilization of the algorithms against real 
data will determine the extent to which the bandwidth must be 
opened to achieve the best results. 

Selection of Q is a problem which perpetually plagues 
users of the Kalman filter. Any given data set will have a 
given Q which will yield best results in terms of a specified 
performance criteria, normally tracking precision. The Q can- 
not be selected with the benefit of hindsight, and must be 
chosen to yield the optimum performance over the ensemble of 
state’ trajectories of interest. The normal method for finding 
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this value would be to take as many raw data tracks as possi- 
ble, and try to determine the Q using statistical methods. 

An alternate technique to generate Q on-line through the 
use of prediction error is described in [10]. In this paper, 
Aldrick and Krabill propose the calculation of Q by the fol- 
lowing method: 

Q(k) = a [Del (k) Del(k) T ] + b[Del(k-l) Del(k-1) T ] (62) 

where a and b are constants to be determined by analysis of 
actual data, and 

Del (k) = X (k/k) - X(k/k-l) . (63) 

This method was investigated to a limited extent, and showed 
some promise if refined. Simple use of (62) seems to open 
the bandwidth wider than is desired. It lias the advantage 
over the use of some constant Q for all runs that, in theory, 
if no uncertainties exist, the Del(k) arrays will be zero 
and thus Q will be 0. Thus, wide bandwidth is achieved only 
when required, as determined by prediction success. In prac- 
tice however, it was found that this method caused the gains 
to oscillate, and created excessive error. 

H. REDUCTION OF COMPUTATIONAL TIME 

The most obvious disadvantage of using the Kalman filter 
compared with the constant gain alpha-beta filter is the in- 
crease in complexity and computation time required to compute 
this optimal gain schedule. There are several techniques 
which can be employed to reduce this burden, all of which re- 
sult in further loss of optimality, but to limited extents. 



35 



i 

In general, the R matrix is a function of the states or 
position of the aircraft. If the rough start and end 

l 

points for the tracks to be followed are known in advance of 
the mission, it is possible to compute the average R, 0, and 
<j> and thus an average R matrix. Unless R, 0, and <p vary over 
wide ranges, the loss in optimality by this approximation 
should be small. Since the gain schedule in state dependent 
only to the extent that R is state dependent, the gain sched- 
ule can now be computed and stored is a mission data table, 

I 

and thus need not be computed on-line. Since 18 different 
i gains are required by RADAR6 and 27 by RADAR9 for each sam- 
pling time, this becomes a rather large problem if auxiliary 
storage is not available. 

I An alternative to storing the gain schedule is to fit 

each of the gain schedules to either an exponential curve or 
to a function of the type 

f (x) =a+bx^+cx^+dx^+ ••• (64) 

Then only 18 (or 27) computations would be required at each 
sampling time. 

Whatever technique is chosen to approximate the true 
gain schedule, it will not be nearly as sub -optimal as a 
constant gain over all time, or two constants which are 
switched in and out. 
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III. PRECISION GUIDANCE SIMULATION 



A. INTRODUCTION 

A limited amount of Precision Guidance simulation program 
documentation was provided with the original version of the 
program. This greatly facilitated understanding of the 
original version of the program, and also subsequent modifi- 
cation for improvement. Many of the original concepts used 
were not changed. Input and output formats were revised for 
ease of use and to provide a greater indication of program 
performance. Since large portions of the original program's 
logic and coding have been retained, these may be mentioned 
and summarized lightly for continuity. The derivations and 
justification of assumptions in these areas will not be ad- 
dressed, as they are considered to be adequately documented 
in the existing system documentation. 

A block diagram of the simulation program control loop 
is presented in Fig. 4. In a few words, the program starts 
with the aircraft at some designated position and velocity 
with respect to a target which is to be bombed. After a 6 
sec filter settling period, control logic determines the 
lateral error which will be incurred if the present flight 
profile is continued until time to release the bomb. Func- 
tions of the lateral error become the controlling mechanism 
to generate command bank angles which cause the aircraft to 
change course. It is emphasized that all control commands 
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Figure 4. Simulation Program block diagram. 



are based on estimated positions and velocities. Thus it 
may be possible to drive the estimated lateral error to zero, 
and still miss the target considerably, if the state estima- 
tion vector is in error. One further comment is that the 
problem is not totally one of estimation and control. The 
aircraft has a bank angle limitation imposed of ± 30 degrees. 
This limits the turning rate such that, depending on the 
initial position and velocity with respect to the target, 
the aircraft may not be physically able to come about to the 
correct heading in the time required. Examples of this situ- 
ation are provided. 

B. INITIALIZATION AND FILTER SETTLING TIME 

The new simulation program is functionally similar to 
the original version in initialization. The same variables 
are used wherever meaningful. The coordinate transformation 
vectors and matrices defined during initialization are re- 
tained, as are all equations for operating in the dive bomb- 
ing mode. (The original program provided logic for executing 
a diving mode, but use of this mode has not been investigated 
in this study.) ‘ _ 

The precision radar sampling rate is 8 Hz. No action ex- 
cept state estimation is taken during the first 2 sec of the 
simulation run. At the 2 sec point, logic is executed to en- 
able prediction of the lateral error at 6 seconds into the 
run. The 4 second lag is assumed to be due to an attempt to 
simulate the fact that the actual computer used is compute- 
time limited, thus constraining the integration logic to 
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I 

execution at no more than 0.25 Hz. At t = 6 seconds, the 
lateral error is first estimated, and the first non-zero 

I 

command to the aircraft can be transmitted. 

C. PRECISION GUIDANCE COORDINATE CONVENTIONS AND TRANSFORMA- 
TIONS 

Three primary coordinate systems are used in Precision 
Guidance calculations. These systems have their origins at 
the target, the radar, and the aircraft. A convention in 
notation is used throughout the program to denote vectors in 

I 

the various coordinate systems. XI represents a position 
vector, XD1 represents a velocity vector, and XDD1 represents 
an acceleration vector. (The only time acceleration vectors 
appear are in RADAR9 . ) The number "1" in the notation simply 
I refers to a specific coordinate system. 

XI represents the state estimate of the aircraft in the 
radar reference frame; X2 represents the true state vector 
of the aircraft in the radar frame. X3 represents the state 
vector of the aircraft in the target frame. Conversion from 
one frame to the other is accomplished through the use of 
transformation matrices EMI, EM2 , and EV1, and the subroutines 
MATMLT , and MATMAD which perform the matrix multiplication 
and addition. The relationship between X2 (or XI) and X3 is 
illustrated in Fig. 5. 

Aircraft control is derived through the use of the X6 
coordinate system, which has its origin at the aircraft and 
its y axis oriented along the estimated aircraft ground head- 
ing. This transformation is accomplished through the coordi- 
nate transformation matrix EM3. Figure 6 .illustrates the 
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North Pole 




The local vertical is labeled Z. Local 
North is the y axis, and local East is 
the x axis in each frame. The dashed 
line is a great circle through the two 
points. Point labeled A/C is the aircraft. 



Figure 5. Illustration of the relationship between the X3 
and the X2 coordinate systems. 
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Figure 6. Illustration of relationship between X6 and the Target. 



the relationship between this coordinate system and the tar- 
get. Note that in the case of perfect estimation and no 
wind that this y axis orientation would point directly at the 
target . 

Each of the X, XD, and XDD vectors are of dimension 3, 
one storage allocation for each physical dimension. Thus, 

Xl(l) = xl 
XI (2) = yl 
XI (3) = zl. 

The coordinate systems and transformation equations are 
fully specified in [11]. 

D. LOOP GEOMETRY AND ERROR CALCULATIONS 

The normal flow of control through the main program loop 
following the initial 6 second settling period begins with a 
time update and a movement of the aircraft in accordance with 
the command bank angle generated at the previous time, all 
within the original subroutine ARCRFT. The true position is 
then transformed into the radar frame, noise added, and the 
new state vector estimated using the RADAR9 subroutine. The 
estimated state vector is then transformed into the X6,XD6 
system, sometimes referred to as the "double primed" reference 
system in the original program documentation. 

As mentioned previously, every 4 seconds the simulation 
program computes new coefficients which are used to calculate 
the ballistic range to the target, and bomb time of fall. 

This involves integration of a system of 16 differential equa- 
tions using a fourth order Rungc-Kutta scheme. The ballistic 
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range (RA) and time of fail for the bomb (IF) are computed 
using a first order linearized approximation to the system 

I 

of equations described above at those times vhen no inte- 
gration has been performed. Subroutines STIFF, DER, and 
OUT are used to perform the required calculations; these are 
completely specified in [10]. 

The ballistic range, RA, is the distance the bomb will 
impact from the present aircraft position. Time of fall, 

TF, is the time in seconds which will elapse between bomb 
release and impact. These values and the ballistic wind 
I components in track and cross-track to the aircraft estimated 
heading determine the impact point for the bomb in X6 refer- 
ence frame. The equations for this estimated impact point, 

I (XGC,YGC), are as described in [11]. Figure 7 illustrates 
the geometry of the situation for the windless case. It is 
assumed that by selecting the bomb release time carefully, 
a near-zero error in impact on the y6 axis can be achieved. 
Only by having the aircraft heading precisely correct will 
the cross track impact error, called lateral error (XE) , be 
zero as well. 

At each sampling point, the Time- to-Go- to-Release the 
bomb, TG, is computed using 

TG = (y 6 - YGC)/y6 . (65) 

The lateral error is given by 

XE = x6 - XGC . 

Other quantities which are used in error calculations are 
the ballistic range to the target, RET, and the heading angle 
error, H1)E. These are given by 
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Figure 7. Illustration of error and bombing geometry. 



RBT = (x5 2 + y 5 2 ) ?s 



(67) 



where X5 is the coordinate system X6 prior to being rotated 
for heading -target alignment, and 

HDE = arcsin (XE/RBT) . (68) 

E. AIRCRAFT CONTROLLER DESIGN 

The new controller for the aircraft is considerably less 
complex than that originally used. The original version em- 
ployed lead-lag networks, suitably digitized, with constants 
which were switched in or out at different stages of the 
simulation run. The lateral error was driven to zero in the 
original version of the program, by selection of a gain con- 
stant times the lateral error to yield a desired bank angle. 

The new controller design attempts to drive both an 
error and an error rate to zero. The error signals to be 
driven to zero are the heading angle error, HDE, and the 
heading angle error rate, HDEDOT, where 

HDEDOT = HDE. (69) 

The desired control bank angle, 4 ) ( j> is given by 

*dl“- G l HDE (70) 

4>d2 = C 2 HDEDOT . (71) 

and 4*^ — 4*^2 + 4*^ 2 ' (72) 

This procedure requires selection of the feedback gain con- 
stants and It should be noted at this time that G ^ 

and G 2 are the only two constants in the simulation program 
which must be determined "through simulation," i.c., by 
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trial and error. The original program was filled with 
numerous "gain constants" which were, or were to have been 

I 

selected by "simulation." In addition, in the original ver- 
sion, much of the theory was developed through the assumption 
that the process was approximately a linear one, and develop- 
ment of a linearized model which reflected these lineariza- 
tions. Such linearization may be somewhat valid in the final 
stages of a long run where only small commands are being sent, 
but in the initial phase of heading correction, the commands 
saturate and the assumption of linearity is invalid. It is 
I in this nonlinear region of operation that commands must be 
optimized to yield a combination of rapid correction of head- 
ing error and minimum overshoot of the correct heading. 

Two sets of control gains have been selected for the 
two prevalent aircraft roll response time constants, t^ = 2 
and t k = 3.3. The gain constants are given below. 

x b = 2 : G ] _ = G 2 = 75 

= 3.3 = G 2 1 : 150. 

These parameters were selected through analysis of numerous 
simulation runs with each and subjectively evaluating the 
resultant performance. The performance criteria used was to 
generate initial commands which are of the "bang-bang" type 
referred to in optimal control theory, causing the bank com- 
mand to saturate at the autopilot limit of ± 30 degrees. 

This causes the aircraft to begin turning toward the correct 
heading at the maximum possible rate permitted. The effect 
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of heading error rate control is to act as a damper in that 
the faster the heading angle changes, the more the control 
will begin to decrease. The second part of the performance 
criteria was to require that when small angle errors are 
noted, only small commands are generated. 

There is a requirement for no bank commands during the 
final second of flight. This is to ensure that the aircraft 
does not make a violent maneuver just as the bomb is dropped, 
causing the bomb to be pitched from the aircraft rather than 
simply dropped as had been assumed in the bonbing calculations. 

For a one second period beginning three seconds before 
bomb release, i.e., TG = 3 to TG = 2, the average estimated 
aircraft bank angle is computed. This average estimated bank 
angle is then sent as the command during the final two 
seconds before bomb release. The reason for sending this 
value vice a command of zero degrees is the possible pres- 
ence of a bias angle. Figure 8 is a plot of a typical simula- 
tion run in which the autopilot has a bias bank angle of 3 
degrees. Note that this bias is automatically compensated 
for by the commands sent. The average command at the end 
of the run is approximately -3 degrees. This has the effect 
of rolling the aircraft back to a level flight profile in a 
steady-state situation, vice a slowly turning profile which 
would be present if the command was 0 degrees. By averaging 
the estimated bank angle over the period defined above, and 
sending this value for the final few -seconds, the aircraft 
continues on an approximately error free path. A zero command 
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at this time would cause the aircraft to begin a slow turn 
and greatly increase the lateral error with no opportunity 
for correction. 

Command system requirements require that the command 
angle be quantized to the nearest 15/128 degree. This is 
accomplished through simple logic as described in [11] . 

F. PRECISION GUIDANCE PROGRAM IMPLEMENTATION 

1. Main Program and Subroutines 

The Precision Guidance simulation program consists 
of a main routine and seven primary subroutines. The pro- 
gram listing is appended at the end of this report. Each of 
the subroutines is included and in each of the routines the 
primary variables and purpose of the routine is specified, 
with the exception of subroutines STIFF, DER, and OUT, which 
serve only to support the Runge-Kutta integration defined 
above. The simulation program was written to run on an 
IBM 360-67 computer system and takes advantage of several 
of the system software subroutines. The program with link- 
age and subroutines requires approximately 150 K to execute 
and will run a typical time of from one to three minutes of 
CPU time. Compilation requires approximately 50 seconds. 

A listing of the variables used in the program along 
with the variable definitions is provided in Appendix A. 

2 . Program Input and Output 

Specific format definitions for program input data 
are provided in Appendix B. The program output has been 
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changed significantly from that in the original program ver- 
sion. A sample of one complete simulation run is provided 
at the end of this report. The initial conditions and con- 
stants are output along with short word definitions of their 
meaning or useage. Each second, a summary of the critical 
parameters of the program are output in block form; the out- 
put key is printed on each page for ready availability to 
the user. During the final second before bomb release, the 
critical parameters summary is printed at every sampling 
point . 

On the final time through the main processing loop, 
all estimates are replaced by their corresponding true values, 
and the bomb "released." The impact point is computed, and 
the miss distances in x,y and overall are printed, designated 
XI, YI , and RI , respectively. 

As a measure of filter effectiveness in noise reduc- 
tion and state prediction in position including deterministic 
motion, the square of filter residual is summed throughout 
the run. The RMS values of filter residue are printed fol- 
lowing the bomb impact miss distances. Filter residue in 
this sense is defined as 

FILRES x = xl - x2 (73a) 

FILRES = yl - y2 (73b) 

FILRES = zl - z2. (73c) 

z 

In addition to the filter residual in each coordinate, 
radial residual defined as 
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a 



FILRES d = (FILRES 2 + FILRES 2 + F1LRES 2 )^ (74) 

K x y z 

is printed in RMS form, where the average is over the entire 

simulation run. 

Some of the critical parameters are stored in arrays 
each second throughout the run. These are then line printer 
plotted at the completion of the run with appropriate 
labelling . 
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IV. COARSE GUIDANCE SIMULATION 

i 

A. INTRODUCTION 

Two documents, [8] and [9], were provided to describe 
the techniques and general approach being followed on the 
original Coarse Guidance simulation program. However, no 
program documentation such as that provided for Precision 
Guidance, [11], was available. The program provided seemed 
inordinately complex in some places and the general approach 

to the problem did not appear to be a viable one from which 

i 

to build an improved version. 

! 

It was determined that the best approach would be to write 
a totally new main simulation driving routine, use the already 

| 

existing aircraft simulation program, ARCRFT, and the already 
developed Kalman filter radar simulation program RADAR6 to 
simulate the Coarse Guidance tracking system. RADAR6 was 
chosen vice RADAR9 simply because neither [8] nor [9] men- 
tioned any difficulty with bank angle biasing. Also, the 
smaller size and faster running time of RADAR6 made that sub- 
routine a preferred choice. The decision to mate RADAR6 to 
Coarse Guidance rather than RADAR9 is by no means a final 
choice. A matter of only a few minutes would be required to 
modify RADAR9 to be compatible with the Coarse Guidance sim- 
ulation program. Thus if biasing is a problem in this mode 
also, RADAR9 could serve as the appropriate unbiased state 
estimator. The changes in RADAR6 from RADAR9 other than 
those which relate to acceleration estimation will be discussed 
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I 

below. A few simplifying changes in the ARCRFT subroutine 
were also accomplished to reduce core storage requirements 

I 

and execution time; these will also be described. 

Some of the concepts and a few equations from the origi- 
nal documents on Coarse Guidance were used in this study. 

Since the new simulation program differs significantly from 
the original version, nearly all equations will be presented 
and most will be derived. 

The basic concept in Coarse Guidance is to get the air- 
craft from some initial starting point to the final bombing 
I run by flying a predetermined course which is specified by 
"waypoints" and azimuths of course "legs." A typical simu- 
lation setup might appear as illustrated in Fig. 9. The 

I 

! designated legs presumably follow a "safe" path for the 
strike aircraft. Also, presumably, if the aircraft deviate 
from the specified path too far, they become in danger. 
Therefore, it is desirable to provide some control to keep 
the aircraft as near to the specified path as possible. Of 
particular importance is to recognize the approach of the 
beginning of a new leg and begin a "command turn" onto this 
new leg at such a time that upon completion of the turn, the 
aircraft will be on the new leg with the same ground heading 
as the leg's azimuth. 

B. PRE-MISSION DATA TABLE COMPUTATIONS AND INITIAL CONDITIONS 
Once the path to be flown is specified, the individual 
legs can be characterized by their azimuth with respect to 
North and their length. The beginning of the mission is 
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Illustration of a typical 3 leg course. 
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Figure 9. 



referred to as a TACAN Entry Point, and is the approximate 
position of aircraft entry into the probles. Wind causes 
the ground heading and the air heading to differ. Since 
most of the aircraft control calculations are made with re- 
spect to the air heading, effects of wind Eust be considered. 
The initial position of the aircraft is placed at the begin- 
ning of the first leg perturbed by some error, and with 
ground heading of the first leg's azimuth also perturbed by 
some angle deviation from desired. 

1 . True and Estimated Wind Components 



wind in the problem. Provision is made for an error in es- 
timated wind speed and direction. All true aircraft motion 
is computed using true wind. All estimated aircraft motion 
and control decisions are made using the estimated wind 
components. Wind is assumed zero in the vertical direction. 
If the true (estimated) direction toward which the wind is 
blowing is 0^ (0 , ) , and the true (estimated) wind speed is 
V (V^) then the components of wind are given by 



Figure 10 illustrates the relationship assumed for 




(75a) 



w, r = V cos(0 ) 
y w v w J 



(75b) 



for the true wind, and 

w xh = V wh sin(e W 



(76a) 




(76b) 



for the estimated wind. 
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North 




Figure 10. Illustration of Wind Relationships. 



North 




Illustration of Waypoint Coordinate Calculations. 
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Figure 11. 



2. Mission Data Table Calculations 



Included in the mission data table are the waypoint 
coordinates, average radar range to each leg, average radar 
azimuth to a given leg, desired air heading while on each 
leg, desired ground speed while on each leg, approximate 
time to fly each leg, and the ground velocity components 
for each leg. 

For notational purposes, it is assnaed that there 
are a total of n legs to be flown; a subscript i on any 
parameter indicates that parameter for the leg. Figure 
11 illustrates calculation of the i+1 waypoint coordinates 
from the previous leg's parameters. If (xwp^ ,ywp^) are the 
coordinates of the i ^ waypoint, the azimuth of the i ^ leg 
is (K , and the length of the i^ n leg is L^, then 

xw P i+1 = xwp i + L i sin(0 i ) (77a) 

y w P i+1 = ywp A + L i cos(0 i ). (77b) 

The average range and azimuth to the leg, and 

Az^ are given by 



R i - h[ (xwp i +xwp i + 1 ) 2 + (ywp i +ywp i+1 ) 2 ]' 2 (78) 



Az . 
i 



arctan 



f xw Pi + xwp. +1 
Lyw Pi + ywp. +1 



(79) 



Figure 12 presents geometry which aids in clarifying 
calculation of the air heading and ground speed, assuming 
the air speed is a known constant. (This assumption is 
maintained through the problem and seems reasonable, since 
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I 




Figure 12. Illustration of Ground Speed and Air Heading 
Calculation Geometry. 
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the pilot has a direct readout of his speed with respect to 
the air.) Since it is desired to fly along the leg, it is 
correct to sum components perpendicular to the leg and set 



these to zero. If V ^ is the estimated air speed and 1L is 
the air heading then 




Summing components in the direction of the leg yields the 
desired ground speed 



The air heading required to fly in the direction 0^ is found 
by solving (80) for IL . 



3 . Initial Position and Velocity of the Aircraft 

The initial true position of the aircraft is that of 
the first waypoint plus a perturbative error. Since the way- 
point is on the ground, the altitude of the aircraft is given 
by the perturbation in the z coordinate. The true position 
of the aircraft is contained in the X array. 




can be broken into Cartesian components, V and as 

follows : 




(82a) 




(82b) 



H. =0. - arcsin 



i i 




(83) 



x3 = xwp. + 

y3 = ywp JL + A 

z3 = A . 
z 



(84a) 



(84b) 



(84c) 
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The initial true ground heading, H , is read as data. 
This in addition to the known airspeed, V , specifies the 

3 . 

true ground speed, V , and the true air heading, H , through 

g a 

a set- of equations similar to (80) , (81) , and (83) . The re- 
sults are given by 



H a = H g " arcsin [v^ sin ( 0 w ' H g )J (85) 

and by 

V = V cos (H - H ) + V cos (6 - H ) . (86) 

g a v a g J w ^ w g 

'The initial true velocity of the aircraft is then broken 

into Cartesian coordinates and stored in the XD3 array. 



x3 = V sin (H ) 
g g 

y3 = V g cos (H g ) 
z3 = 0. 



(87a) 

(87b) 

(87c) 



Note that it is through specification of H different from 

g 

6., and A and A different from 0 that nonzero perturbations 
1 ’ x y r 

in velocity and position from the desired values are entered. 
Figure 9 also shows the geometry which might exist in the 
case of nonzero displacement and velocity from the desired 
track. 



C. AIRCRAFT POSITION, VELOCITY, AND ERROR ESTIMATION 

As previously stated, aircraft motion is simulated by the 
use of a slightly modified version of subroutine ARCRFT as- 
sociated with Precision Guidance. Position and velocity es- 
timation is accomplished using RADAR6 . At the beginning of 
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the main processing loop, true aircraft position is updated 
followed immediately by a prediction update on estimated 
aircraft position. A new estimation update is performed 
only if the total radar sampling interval DTRAD has ellapsed. 
(It is assumed that prediction updates occur at a higher 
rate than radar sampling.) The most current aircraft posi- 
tion and velocity estimates are contained in the XI and XD1 
arrays. From this and the estimated wind components, the 
estimated air and ground headings, and estimated ground speed 



are computed; H 


H , , and 
gh 


V gh ’ 


respectively . 




H , 
gh 


= arctan 


M 




(88) 


H , 


= arctan 


"xl - 


W xh 1 


(89) 


ah 


L 


-yi - 


Wy h J 




V u 
gh 


= ( xl 2 + 


yi 2 ) % 


* 


(90) 



These parameters can be monitored to determine the degree 
of error in heading and speed from the desired values. An 
additional and important parameter to be monitored is the 
extent to which the aircraft has deviated from the desired 
course, E ^ . Figure 13 illustrates the geometry used in 
this calculation. The error is given by 

E es t = ( xl " xw Pi) cos(6 i ) - (yl - ywp i ) sin(O i ). (91) 

Note that (91) computes the distance from the present posi- 
tion estimate to the leg i. In a command turn, the aircraft 
should not be on either leg, but wilJ be somewhere between 

0 2 




Figure 15. Illustration of Geometry for Calculating Deviation 
from t lie Desired Leg Path. 



6 3 



til 

the legs. For this reason, distance from both the 1 and 
i+1 leg is computed, and the smaller of the two values 
chosen . 

True headings, speeds, and distances corresponding to 
those given above are also computed using the same equations 
with X3 substituted for XI, and W substituted for W^. 

D . COMMAND TURN CALCULATIONS 

The command turn is that which is computed to cause a 
smooth transition of the aircraft from the present leg to 
the next leg with a minimum or overshoot or undershoot of 
the desired path. The aircraft autopilot is constrained to 
a maximum bank angle, typically equal to 30 degrees. 

Since the heading angle rate is proportional to the bank 
angle in a coordinated turn, the time to complete the turn 
using the maximum possible bank must be precalculated so that 
the turn can be started prior to reaching the new leg. This 
process is complicated by three factors. The first is that 
the time to the next leg is variable with the current posi- 
tion and velocity. Thus the time to begin the turn as well 
as the amount of heading change required is a function of the 
state estimate. The second factor is that the equation to be 
solving for the amount of time required for the turn is a 
transendcntal equation and must be solved through iteration. 
The third factor is that this equation becomes unduly complex 
if the turn is not started from a zero bank angle. This rc- 
qui res that the aircraft begin the turn in level flight. 
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1 . Time Remaining on Present Leg 



The time remaining on the. present leg before inter- 
secting with either the next leg or its extension is com- 
puted by finding the intersection of the two paths, calculating 
the distance to be traversed, and then dividing by the esti- 
mated ground speed. The geometry is illustrated in Fig. 14. 

The coordinates of intersection of the present path, 
based on the present velocity estimate and the next leg are 

designated (x^ >Fi n t^ * Let m i anc * m i + l re P resent the slopes 
of the present aircraft heading and the next leg, respectively. 
Then 

itk = yl/xl (92) 

= (ywp^i - ywp.;)/(xwp - xwp.). (93) 

_L • -1 J. ■ -L _L -k. ■ -L -L 

The equations of the two lines whose intersection is to be 
found are 

y = m i+1 x - (m i+1 xwp i+1 - ywp. +1 ) (94) 

y = m^ x - (m^ xl - yl) . (95) 

These equations are solved simultaneously to give the point 
of intersection. 



x- . 
mt 



(yl - ywp i+1 ) + (m i+1 xwp i+1 - m i xl) 
(m i + i - mfj 



(96) 



y. . = m. x. . 
'int i int 



The distance between xl 



m. xl + yl. 
i 1 

and the point of intersection 



D = f (x . - xl) 

tg 1 v ml J 



Ojnt ‘ 



( 97 ) 

is then 
(98) 
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Figure 14. Illustration of Path-Leg Intersection Geometry 



i + 2 
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and the time to reach this point is given by 

TLEG1 = D /V , . (99) 

tg gh 

2 . Time Required to Complete the Turn 

The amount of turn required, AH, is simply 



AH = |H. +1 - H ah i . (100 ) 

Air headings are used in the calculations since all aircraft 
motion and turn equations must account for the possible pres- 
ence of wind. Equation (50) gives the relationship between 
the change in heading angle, time in the turn, and the roll 
response of the aircraft. Setting A^ equal to AH and 4> 

( 101 ) 



equal to <j> gives 



r. 



AH = (g/V ah ) i <p m T + (ii>(k-l) - 4» m ) (T b ) (: 



; T/ xi 



If it is assumed that <+> (1c - 1 ) is zero (starting with level 
flight) , then the equation can be rewritten as 

K = U - (1 - e' U ) (102) 



where 



U = T/t, 



(103) 



and 



K = 



AH V , 
ah 

g T, O 

6 b m 



(104) 



For the purpose of finding the turn time required, 
assume that the turn can be divided into three parts as il- 
lustrated in Fig. 15. The first part is the transient period 
during which the aircraft is coming to the maximum (or minimum) 
bank angle, going to that limit exponentially with time constant 




Figure 15. Illustration of Relationships Between Times and 
Bank Angles in a Command Turn. 
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The second part of the turn is that segment when the 
aircraft is holding at the bank limit and changing direction 
with an approximately constant heading rate. The third and 
final part of the turn is a transient segment during which 
the aircraft returns to a zero bank angle, also exponentially. 

Figure 15 shows that the turn can be divided into 
two approximately equal parts, during which the aircraft 
executes approximately half of the turn, AH/2. Let the time 
required to complete half of the turn be T ^ > an< 3 the time to 
complete half of the constant bank segment of the turn be 
T^. If K, above, is replaced by K/2, then the solution to 
(102) yields 

u - T t2 / V C105) 

The equation is solved using the Newton-Raphson iterative 
technique. The initial approximation to T 2 is taken from 
[9] . Iteration continues until the change in T 2 is less 
than 0.01. 

Motion of the aircraft through the turn is approximated 
in nearly the same manner as described in [9], Briefly, it is 
assumed that the aircraft continues on an approximately a 
straight path for a period ^+2 ^02 secon< ^ s after the command 
turn is ordered, followed by the command turn as described 
above, followed by another period of approximately straight 
flight along the new leg during the final seconds of the turn. 

A diagram of the geometry of the turn is presented in Fig. 16. 
Since the turn has been computed with respect to the motion 
of the air mass, movement of the air mass during the turn is 
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Figure 16 . Illustration of Turn Geometry, Including Effects of Wind Shifting the Final Aircraft 
Position. 



compensated for. The only modified equation used to calculate 
the distance prior to leg intersection from which the turn 
must start is that which relates to air mass motion. The new 
equation for computation of the parameter D c is given cor- 
rectly by 

D = T , V , s in (6 . - - 9 ,) sin(6. , - 9.). (106) 

c o2 wh v i+I wlr v i+I i J ^ 

The turn must be started a distance D, before the intersec- 

o 

tion with the next leg, where is shown on the illustration 
and mathematically defined adequately in [9] . 

3 . Command Turn Timing Logic 

The time-to-go before beginning a command turn is 

given by 



TG = TI.EG1 - IK/V , 
o gh 



(107) 



When TG is less than or equal to zero, the aircraft is com- 
manded to go to the maximum bank angle, with the sign of the 
bank chosen appropriately. The time to command the bank 



angle back to zero is T st0 p tn anc ^ I s given by 

T . . = T - 2t, . 

stoptn t b 



(108) 



Initially it might appear that 3r^ should be subtracted from 
the total turn time, since that would be closer to the amount 
of time required to decrease the full bank angle when the 
change is occurring exponentially. However, this value was 
arrived at through simulation trials, and is probably best 
due to the compounded approximations made in the overall turn 
solution. The time at which the turn is complete is defined 
as T ; no actions depend on this time. 
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Figure 17 is a logic flow chart of the start/stop 
turning process. Briefly, a counter in the form of the 
variable TINTRN is incremented each time through the loop 
as the turn progresses. The variable ITURN is used as a 
flag to pass logical control to the correct coding. When no 
turn is in execution, ITURN = -1, and the only turning which 
is performed is that to correct the course deviations. No 
control banking for the above purpose is permitted within 
the 3 t^ seconds prior to executing a command turn. This is 
to ensure that the bank angle of the aircraft is zero when 
beginning the command turn, an assumption which was used in 
the derivation of the command turn equations. During the 
command turn, ITURN = 0, and all control logic is bypassed. 
Immediately upon executing the logic which indicates the turn 
is ending, ITURN = 1, and calculations to check for the next 
turn time begin. If required, a new command turn can be 
executed immediately, before the aircraft has come down from 
its bank from the previous turn. This feature was added to 
ensure that if a short leg or a very acute turn was encoun- 
tered, the best possible flight trajectory would be flown. 
Examples of the requirement for this feature and its perform- 
ance are included. 

E. AIRCRAFT CONTROLLER DESIGN 

The controller, used to guide the aircraft to a desired 
velocity and keep it there until a command turn, is of the 
same type used in Precision Guidance. In this case, the con- 
troller attempts to keep the estimated ground heading, II 
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Figure 17. Command Turn Start and Stop Logic. 



at that heading from the current estimated position which 
will fly the aircraft directly at the next waypoint. This 
scheme was selected for two primary reasons. The first is 
that it requires minimal control to get the heading correct 
and thus requires considerably less updating than the tech- 
nique used in the original simulation program. The second 
reason is that this technique tends to avoid the problem of 
oscillation which plagued the original scheme. The control- 
ler features a technique which will block any command updates 
to the aircraft unless the heading error exceeds some minimum 
error angle, H . . 

The estimated heading to the next waypoint is 



H, . — arc tan 
hgwpt 



r ~Pi»i 

L>" w Pi.i 



xl 

Ti_ 



(109) 



The heading error is then 



HUE = H , - H, 

ghwpt hg 



( 110 ) 



and the heading error rate, HDF.DOT is the discrete derivative 
of (110) . These are combined as given in (70) , (71) , and 

(72) to give the desired control angle. The gains G-^ and G ? 
were selected equal,, as before, with the value 4. Different 
feedback gains were used in a number of simulation runs, and 
the results indicated that the system was somewhat insensi- 
tive to the values chosen. The gains yielded desired commands 
of about 20 degrees during normal course corrections, and with 

the H . feature seemed to avoid the undcsired oscillation 
crmin 

for the most part. 
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The commands are quantized before being sent, as in 
Precision Guidance. 

F. COARSE GUIDANCE PROGRAM IMPLEMENTATION 
1 . Main Program and Subroutines 

The Coarse Guidance simulation program consists of 
a main routine and two primary subroutines. The subroutine 
ARCRFT used to simulate true aircraft motion, was simplified 
from that version used in Precision Guidance in that all 
dive bombing equations were removed. The RADAR6 subroutine 
differs from RADAR9 by the obvious fact that it does not es- 
timate acceleration due to bank angle bias, as well as in two 
other ways. The first is that RADAR6 is designed to operate 
at a sampling interval greater than the control interval. 

This requires that the subroutine be called, standard predic- 
tion equations executed, and then estimation only if the sam- 
pling interval has elapsed. The second difference is in the 
use of NWLD. NWLD indicates that "wild points" are in ef- 
fect, and are ignored by the radar filter. This amounts to 
simple prediction without the benefit of radar sampled data. 

The Coarse Guidance radar receives little if any re- 
turn from the aircraft during command turns. In an effort 
to simulate this effect, NWLD is set to a negative value dur- 
ing turns, thus requiring prediction of position and velocity 
based only on past data and commands. 

Coarse Guidance requires a little less than 100 K of 
core for linkage and execution. Compilation of the program 
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and subroutines requires about 25 seconds, and execution 
requires between 5 and 25 seconds, depending on the number 
of legs and the leg lengths. 

A listing of the variables used in the program along 
with the variable definitions is provided at the end of 
Appendix A. 

2 . Program Input and Output 

Specific format definitions for program input data 
are provided in Appendix C. 

The program output consists of three primary parts. 
The first part is a listing of all input data and simulation 
run parameters. The second part is a summary of critical 
parameters printed once per second, and the third part is a 
summary of the deviation error from the desired path, along 
with a plot showing the desired, true and estimated paths. 

At the beginning of each new leg, a summary of data 
table information is provided about that leg which can then 
be compared to what actually occurs. As the run progresses, 
a notification of the three stages of command turning is 
presented interspersed with the critical parameter summary. 

The only performance criteria selected is the root 
mean square aircraft deviation from the desired path. The - 
plot of the paths shows the deviations very well. The plot 
scaling is accomplished to maintain the same scale on both 
axes for a truer visual indication of relative error. All 
output values are in nautical miles and feet per second, for 
position and velocity, respectively. 
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V. PRESENTATION OF RESULTS 



A. PRECISION GUIDANCE PERFORMANCE COMPARISON 

In the process of improving, changing, and verifying per- 
formance of the new version of the Precision Guidance simula- 
tion program, a very large number of runs were accomplished. 

A proper comparison between the two program versions re- 
quires side-by-side contrast of appropriate performance pa- 
rameters, each derived from the same set of initial conditions. 

In establishing specific run parameters for the initial 
conditions, three parameter sets were varied: initial air- 

craft position, X3, autopilot bias angle, 4>^ , and aircraft 

vq! 1 ro S130 ri c ° i r'onc ^ on f t All nflior 3.t) "I uc i -n p> 

^ ~ ~ ^ • 

programs were held constant and equal in the two program ver- 
sions. Variation of the above parameters seemed to yield 
conditions which would ably show those areas where perform- 
ance was improved, and at the same time permit runs with 
initial conditions varying from nearly perfect to considerably 
in error from the optimal initial bombing path. 

Ten different aircraft initial positions were chosen and 
assigned a "run number." The initial velocity on each of 
these was identical. The resultant geometry created is shown 
in Fig. IS. Note that runs 1 through 5 present a consider- 
ably more difficult mission than runs 6 through 10. "Diffi- 
culty" can be roughly equated to the angle through which the 
aircraft must change its velocity in order to fly toward the 
target, located at the origin of the X3 coordinate system. 
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1 . Filter Performance Comparison 



The filters were tested within the framework of the 
Precision Guidance program. The most important performance 
parameter to monitor for this problem is the difference 
between the estimated position of the aircraft and the true 

position of the aircraft, already defined as the "residual" 

* 

in II and III. Table I presents a comparison between the 
Kalman filter and the alpha-beta filter's radial residual, 
i.e., total estimation error in position, and Fig. 19 pre- 
sents these results graphically. The numbers shown are aver- 
age radial residuals, averaged over the runs accomplished for 
each version. 



AVERAGE RADIAL RESIDUAL (ft) 



Init i al 
Run No. x5 

1. -50,000 

2. -40,000 

3. -30,000 

4. -20,000 

5. -10,000 



Original Version 
79 
107 
154 
216 
128 



Improved Version 
15 
17 
20 
21 
27 



\ 

Table I. Comparison of Average Radial Residuals for the 
Alpha-Beta and Kalman Filters. 



The figures show an average percent improvement for 
the new program version of 680 percent, with a maximum im- 
provement of over 1000 percent. It should be noted at this 
poin.t that these numbers do not represent the results of a 
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SO 



true Monte Carlo simulation. To achieve Monte Carlo pre- 
cision would have required an inordinate amount of computer 
time to prove or illustrate a point. The differences in ac- 
curacy presented in Tablclare not simply a result of stochas- 
tic luck on obtaining a "good" set of random numbers . The 
differences in the two versions represents a biasing in the 
original filter due to the lack of deterministic forcing. 

This point is brought out even more forcfully in the com- 
parison of Figs. 20 and 21, and Figs. 22 and 23. These plots 
show that the errors in estimation for the Kalman filter on 
a typical run are roughly unbiased; hox^ever, the biasing on 
the estimation in both x and y coordinates from the alpha- 
beta filter is very obvious. Note that, as explained in II, 
the alpha-beta filter eventually "catches up" to the correct 
position. This occurs only after the initial period of 
banking at the limit is complete, and explains the peak in 
residual error which can be noted on the runs beginning at 
-20 K ft. 

In theory, even if the alpha-beta filter had effects 
of deterministic forcing included, the Kalman filter would 
perform in a superior manner, due to its having an "optimal" 
gain schedule to give a minimum state covariance of error. 

2 . Bombing Accuracy and Time Response Comparison 

Two performance parameters on the overall Precision 
Guidance routine were observed and minimized throughout 
program development. The importance of actual bombing ac- 
curacy is obvious. A less obvious' but very important parameter 
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is the time required on the final bombing leg in order to 
achieve "sufficient" accuracy. The longer that a pilot is 
required to fly a straight path toward a target, the more 
his chances of being shot down increase. It is then desir- 
able to swing onto the final leg from Coarse Guidance and 
complete the mission as soon as possible. 

This time, to be designated T, was difficult to mea- 
sure directly. It is a strong function of the amount of 
heading angle change required from Precision Guidance ini- 
tialization. By varying the angles and distances to the 
target -in the runs (1 through 10), it is possible to observe 
the amount of time required for the lateral error to reach 
some threshold minimum acceptable value, under which, in an 
absolute value sense, the bomb impact will be scored a suc- 
cessful mission. For the purpose of this study, this value 
was arbitrarily set at 70 feet. 

Table II presents a summary of data from the simula- 
tion runs illustrating bombing accuracy and the time response 
parameter T. Note that most of the lists of bombing ac- 
curacies have a letter designation at the top of the respec- 
tive columns. There are referenced on Figures 24 and 25, 
where comparisons of the bombing accuracies are presented 
graphically . 

The inaccuracies resulting from those runs which start 
relatively close to the target are due to the fact that the 
aircraft cannot turn at a high enough rate to get on a bombing 
line passing near the target before it must drop the bomb. A 
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Precision Guidance Simulation Results. RI is bomb impact miss distance. T is time 
into run at which the lateral error dropped below 70 ft. 



Figure 24. Illustration of Improvement in Bomb Impact Miss 
Distance for Varying Run Lengths. Initial y3 = 
20,000 ft. 

Original Program Results: Curves A,B 

Improved Program Results: Curves C,D 




-10K -20K - 3 OK -40K -50K 

Feet 
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Figure 25 



10000 - 



Illustration of Improvement 
in Bomb Impact Miss Distance 
for- Varying Run Lengths. 
Initial y3=5000 ft. 

No bias angles. 



Bomb Impact Miss Distance 
vs . 

x3 Simulation Run Length 



0 Original Program Results (E) 
Improved Program Results (F) 




-10K 



-20K 



- 4 0 K 
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large part of the error improvement which can be seen in 
these two figures is due to the filter improvement. However, 
the controller also has a large effect on the overall re- 
sponse of the aircraft. The use of angle rate feedback helps 
to drive the lateral error to zero faster than in the origi- 
nal version. This is presented graphically in Fig. 26. Note 
also that while the original version suffers the effect of a 
large overshoot past zero lateral error, the improved version 
does not. This is fundamentally the result of the rate feed- 
back acting as a damper. 

Examination of T in Table II whill show that the or- 
iginal simulation program achieved a lateral error below 70 
feet only once in the simulation trials run. The improved 
version reached the acceptable error 9 out of 15 times. In- 
spection of the times will show that very little time is 
required to achieve this figure of accuracy, including those 
runs which are initialized with rather extreme initial error 
angles (e.g., 3, 8, and 9). Use of the new algorithms to 
estimate position and control the aircraft should result in 
a marked improvement in overall accuracy in bombing and in 
time on the final leg required to achieve this accuracy. 

B. COARSE GUIDANCE PERFORMANCE 

No contrasting performance tests were performed for the 
Coarse Guidance mode. The primary reason for this is the 
lack of sufficient program documentation on the original ver- 
sion of Coarse Guidance. Thus, the Coarse Guidance results 
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LATERAL ERROR (XE) COMPARISON PLOT 
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6,00 26,13 46.25 66.38 06.50 



presented must stand on their own merit, on an absolute 
instead of a relative performance scale. 

As was the case for Precision Guidance, a large number 
of Coarse Guidance simulation runs were completed in the 
process of obtaining the final program version. A sample 
of some of the representative track plots with the RMS track 
deviations for those runs arc presented in Figs. 27 through 
30. The significant aspects of each of these runs are 
discussed below. 

Figure 27 represents a rather simple track to follow, and 
is probably the most realistic track of the four presented. 
The turns to be executed and the initial heading and position 
errors are not excessive. The control and turn logic caused 
the aircraft to follow the desired path very closely. 

The scale on Fig. 28 is about half of that for Fig. 27 
and thus the initial errors are more obvious. The track 
deviation figure is higher here than before due to the fact 
that the initial errors tend to take the aircraft away from 
the desired track, vice in the track's general direction, 
as was the case for Fig. 27. Note that the control logic 
has caused the aircraft to fly toward the second waypoint. 

At the appropriate time, the aircraft begins a command turn 
in order to be aligned properly on the second leg upon turn 
completion. The error from this point on is minimal, never 
exceeding 0.28 NM on the final leg. 

Figure 29 illustrates a most unlikely track to be en- 
countered in practice, but one which checks oul the program 
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DESIRED, TRUE, AND ESTIMATED POSITION 
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algorithms quite well. The initial errors are as before. 

In this case, the aircraft exists from the first command 
turn only to note through turning logic that it is already a 
few seconds behind schedule to be able to make the correct 
turn onto the final leg. This is because of the combined 
condition of a relatively short leg length remaining and a 
requirement to execute 135° change in heading. It can be 
seen from the plot that the aircraft was not able to come 
to the correct final leg heading of 225° from the command 
turn control alone. Another effect which is evident from 
this run is the lack of position update information in the 
filter during a command turn, due to poor radar return 
simulation . 

Figure 30 is the result of a track of the same shape as 
that used in Fig. 29, except the legs are twice as long. 

The initial position error is also the same as that for 
Fig. 29, but the initial heading error from the desired 
track in this case is 90°. The purpose of this track is to 
show that the program's algorithm can handle effectively a 
case such as that in Fig. 29, provided the given physical 
constraints of the system permit time to react. 

Note that the aircraft performs the correct course change 
and heads for the second waypoint, and from that point on 
is nearly coincident with the desired track. The relatively 
high track deviation figure for this run is due almost en- 
tirely to the effects of the initial heading and position 
errors . 
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Each of these runs was made with a minimum error heading 

for correction, H . , of 5°. Use of this number on runs 
* ermm’ 

for which excessive initial errors were not present, i.e., 

Figs. 27 and 28, resulted in a command correction being sent 
less than 10 percent of the time during which command turns 

and initial error corrections were not in progress. This 

* 

indicates a vast improvement towards the goal of decreasing 
the numbers of mid-course corrections being sent to the 
pilot during this phase of the problem. 

Autopilot control of the aircraft during this phase was 
assumed because it was not known how to realistically model 
the effects of a pilot receiving an instruction to come to a 
new heading. If pilot control of the aircraft is desired, 
the algorithms used provide a heading error which could be 
used. In fact, it is this value which is used to compute 
the bank angle correction. 

In summary, the improved Coarse Guidance algorithms ap- 
pear to yield a greatly simplified and accurate technique 
for guiding the aircraft onto the final leg. In addition, 
the errors in position and heading as the aircraft enters 
the final leg are such that the Precision Guidance filter 
can be initialised with valid nonzero velocities, and thus' 
the time required on the final leg can be reduced even further, 
since the requirement for a full six seconds without control 
could be reduced. 
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APPENDIX A: LIST OF PRINCIPLE VARIABLES FOR THE PRECISION 

AND COARSE GUIDANCE SIMULATION PROGRAMS 



A. PRECISION GUIDANCE MAIN ROUTINE 



ATITLE-ZTITLE 

BIS 

BFF 

D 

DEG 

DT 

DTCON 

EMI 

EM 2 

EM3 

EV1 

FILRES 

G 

HDE 

HDEDOT 

HA 

IB1 

I DTCON 

ITAB1 

ITH 

IU 

KTH 

MWLD 



Titles of the plots 
Means of RADAR9 R, 0, and <f> noise 
Bomb ballistics form factor 
Bomb diameter 

Radian/degree conversion factor 

Radar sampling interval 

Control update interval 

Matrix to convert from target to 
radar reference frame 

Matrix to convert from radar to 
target frame 

Matrix to rotate from X5 to X6 
reference frame 

Translation matrix to go from 
target to radar reference frame 

Contains estimation error in 
radar frame for x,y,z, and radial 

Gravitational constant, 32.2 ft/sec 2 

Heading error angle 

Time derivative of HDE 

Height of target above sea level 

Equal to 0 in level bombing mode 

Integer relating DT to DTCON 

Counter for number of points to plot 

Equal to the number of times through 
main processing loop, except on last 
time through, when equal to -1 

Seed for random number generator, 
NORMAL 

Counter to indicate time for 
integration for RA and TF calculations 

Number of "wild points" to be thrown 
out by RADAR9 as invalid 
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NWLD 


Set equal to MWLD at time TWLD and 
decremented as "wild points" are used 


NTB 


Number of bombing constant table to 
be used 


PH 


True bank angle 


PHI 


Estimated bank angle 


PHD 


Desired feedback command 


PHD1 


Desired feedback command from HDE 


PHD2 


Desired feedback command from HDEDOT 


PHDAVG 


Command to be ordered for last two 
seconds before bomb release 


PHB 


Autopilot bias bank angle 


PHC 


Command bank angle sent 


PS 


True aircraft heading angle in 
target frame 


PS1 


Estimated aircraft heading angle in 
radar frame 


PSD 


True aircraft turning rate 


RE 


Radius of the Earth 


RBT 


Distance from aircraft (estimated) 
to the target 


SIG(l) ■ 
SIG (4) - 
SI GW 


- SIG(3) Standard deviations of R, 0, and <J) 

■ SIG(6) Initial radar velocity estimates 

Standard deviations of random forcing 
to be assumed for calculation of Q 


T 


Time into the simulation run 


TB 


n j the aircraft roll response time 
constant 


TG 


Timc-to-Go to release bomb 


TF 


Time of fall for the bomb 


TLVL 


Time for required level flight before 
bomb release; equal to two seconds 


TLVL1 


Equal to TLVL +1 


VE 


Total airspeed 


VEH 


Horizontal airspeed 


W 


Weight of bomb in pounds 


WH 


Estimated wind in target frame 


WR 


Estimated wind in radar frame 


WT 


True wind in target frame 



1.00 



XI , XD1 ,XDD1 


Estimated position, velocity and 
acceleration in radar frame 


X2 ,XD2 


True position and velocity in 
radar frame 


X3,XD3 


True position and velocity in 
target frame 


X5,XD5 


Estimated position and relative 
velocity of target with respect 
to the aircraft; "primed system" 


X6,XD6 


X5,XD5 system rotated to align 
the YD6 axis pointing at the 
target; for printing and plotting 
purposes, this vector contains the 
error between true and estimated 
position in the X3/X5 frame 


XE 


Lateral error 


XGC 


X bomb impact point in X6 frame 


XXA-XXG 


Arrays used to store variables for 
plotting at program end 


YYA-YYU 


Arrays used to store variables for 
plotting at program end 


YGC 


Y bomb impact point in the X6 frame 


B. SUBROUTINE ARCRFT 




CA1 - CA5 


Constants used in aircraft motion 
equations 


DT 


Radar sampling interval; also interval 
of update for aircraft true position 


DEG 


Radian/degree conversion factor 


DT3 


Equal to DT/2 


G 


Gravitational acceleration, 32 ft/sec 2 


ITH 


Equal to number times through loop, 
except on last time when ITH = -1 


IB1 


Equal to zero for level flight mode 


PH 


True bank angle 


PHB 


Bank angle bias 


PHC 


Command bank angle 


PHN 


New bank angle 


PS 


Heading angle 


PSD 


Heading angle rate 


PSDN 


New heading angle rate 
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PSN 


New heading angle 


T 


Elapsed time since start of run 


TB 


, the roll response time constant 


TG 


Time to go to release bomb 


VT 


Horizontal airspeed 


WT 


True wind vector 


X3 


Aircraft position in target frame 


XD3 


Aircraft velocity^ in target frame 


C. SUBROUTINES 


RADAR6 AND RADAR9 


A 


Azimuth angle 


ADDSUB 


Subroutine entry to add or subtract 
two matrices 


ADUM 


Temporary array used in matrix 
arithmetic 


ANGMAX 


Maximum angle, used to prevent 
overflow 


ANGMIN 


Minimum angle, used to prevent 
underflow 


BIS 


Range, azimuth, and elevation noise 
bias 


BDUM 


Temporary array used for matrix 
arithmetic 


CA1 ,CA4 , CAS 
CAA2 


Constants used in aircraft motion 
equations 


DEG 


Radian/degree conversion constant 


DT 


Radar sampling interval in RADAR9 


DTRAD 


Radar sampling interval in RADAR6 


DTCUM 


Accumulator to determine next time 
to sample position in RADAR6 


DT2 


Equal to DT 2 /2 


DTRAD 2 


Equal to DTRAD 2 /2 


DT3 


Equal to DT/2 


DT4 


Equal to DT 3 /6 


E 


Elevation angle 


E 1 , E 2 ,E3 


Difference between predicted position 
and data sample during state estima- 
tion 


G 


Gain matrix 
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GG 

GAMMA 

GNX,GNY, GNZ 

H 

HT 

INVERT 

ITH 

IU 

NORMAL 



PE 



PHI 

PHC 

PHI 

PHITRN 

PROD 

PS1 

PSN 

PSD1 

Q 

R 

RAN 

SIG(l) - SIG (3) 
SIG (4 ) - SIG (6) 

SIGN 

TB 

TRANS 



Gravitational acceleration, 32 ft/sec 2 
F matrix 

G(l,l), G(4,2), and G(7.3) terms o£ 
the gain matrix 

Measurement matrix 

Measurement matrix transposed 

Subroutine entry to invert a matrix 

Equal to number of times through 
main processing loop, except on 
last time through when ITH = -1 

Seed for random number generator, 
NORMAL 

Subroutine to generate normally 
distributed random variables to be 
used as noise 

Covariance of estimation error 
matrix 

Estimated bank angle 

Command bank angle 

State transition matrix used for 
gain generation only 

Transpose of PHI 

Subroutine entry to perform the 
product of two matrices 

Estimated heading angle with respect 
to the wind 

Predicted new heading angle before 
acceleration estimate correction 

Estimated turning rate 

Array resulting from expected random 
forcing assumption 

Range to target 

Contains random numbers from NORMAL 

Standard deviation for R, 0, and <}> 

Initial velocity estimates for 
prediction 

Standard deviations of random forcing 
in x,y, and z of radar frame 

Roll response time constant, 

Subroutine entry to perform the 
transpose of an array 
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VT1 

VARR , VART , VARP 
WR 

X1,XD1 ,XDD1 
X2,XD2 

X1P ,XD1P , XDD1P 

XDATA 

XIDENT 

D. COARSE GUIDANCE 
ATITLE-ETITLE 
AZ 

AZI 

BIS 

DEG 

DEL 

DT 

DTCON 

DTRAD 

DTG 

EEST 

ETRUE 

FEET 

G 

G1 

G2 

H 

HEA 

HTA 

HEG 



Estimate of horizontal airspeed 

Variances of R, 0, and <f> 

Estimated wind vector in the radar 
frame 

Estimated state vector in radar 
frame (XDD1 not used in RADAR6) 

True state position and velocity 
in the radar frame 

Predicted state vector in radar 
frame 

Cartesian radar observation with 
noise 

Identity matrix 

MAIN ROUTINE 

Titles for plots 

Array for storing average leg 
azimuths 

Initial azimuth for first leg 

Array for RADAR6 noise biases 

Degree to radian conversion factor 

Initial displacement vector for true 
aircraft position 

Set equal to DTCON 

Control update interval 

Radar sampling interval 

Distance from XI to leg intersection 

Estimated aircraft deviation from 
leg 

True aircraft deviation from leg 

Feet - to-nautical mile conversion 
factor 

Gravitational acceleration constant 
HDE gain constant, G^ 

HDEDOT gain constant, G^ 

Array of desired air headings 
Estimated aircraft air heading 
True aircraft air heading 
Estimated aircraft ground heading 
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HTG 




True aircraft ground heading 


HDE 




Heading angle error 


HDEDOT 




Heading angle error rate 


HERMIN 




Minimum error angle for which 
commands will be sent 


I 




Equal to the index on the current 
leg 


IDTRAD 




Integer relating DTCON and DTRAD 


IEND 




Equal to 0 except on last leg 


ITH 




Equal to the number of times through 
the main processing loop 


IU 




Seed for the random number generator 


NLEG 




Number of legs for the run 


PHB 




Bias angle; assumed zero 


PHC 




Command bank angle 


PHD 




Desired bank angle from controller 


PLENGTH 




Length of each leg 


RAD 




Radian to degree conversion factor 


RMSEST 




RMS value of EEST 


RMSETR 




RMS value of ETRUE 


RANGE 




Array containing average range to 
each leg 


RANGE I 




Range to start of first leg from the 
radar 


SIG(l) - 


SIG (3) 


Standard deviations for R, 0, and <j) 


SIG (4) - 


S I G (6) 


Initial velocities for filter 
prediction 


SI GW 


- 


Standard deviations for assumed ran- 
dom forcing; assumed zero for all 
runs in this study 


T 




Time into simulation run 


TB 




True aircraft roll response time 
constant 


TAUH 




Estimated aircraft roll response time 
constant; set equal to TB for all 
runs in this study 


TINTRN 




Time into a command turn 


THETA 




Array containing azimuths of each leg 


THETW 




Direction toward which true wind blow: 



105 



THETWH 

TLEG 

TLEG1 

TSTOP 

TSTPTN 

TTOTTN 

TG 

VA 

VAH 

VEG 

VTG 

VW 

VWH 

WR 

WT 

XI ,XD1 

X2 ,XD2 

X3,XD3 

XWP 

XXA-XXE 

YYA-YYE 

YWP 

Note: There are 

variables 
except the 
miles . 



Direction toward which estimated 
wind blows 

Array containing desired time to 
cover each leg 

Estimated time remaining on a leg 

Time at which problem ends 

Time at which a zero bank command 
is sent to end a turn 

Time at which a turn is completely 
over 

Time to go to start a turn 
True airspeed 
Estimated airspeed 
Estimated ground speed 
True ground speed 
True wind speed 
Estimated wind speed 

Equal to WH; components of estimated 
wind 

Components of true wind 

Estimated aircraft position and 
velocity 

Desired track data; used for plotting 
only 

True aircraft position and velocity 

Array containing x waypoint co- 
ordinates 

Arrays for use in plotting 

Arrays for use in plotting 

Array containing y waypoint co- 
ordinates 

several variables with the suffix NM. These 
correspond to the variable without the suffix, 
NM indicates the value stored is in nautical 
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APPENDIX B: PRECISION GUIDANCE INPUT DATA FORMAT DEFINITION 



CARD 


COLUMNS 


FORMAT 


VARIABLE 


REMARKS 


1 


1-10 


F10.3 


WT(1) 


True wind in target frame 




11-20 






WT ( 2 ) 


x,y, and z direct ions ;WT (3) 




21-30 






WT ( 3) 


is usually zero 




31-40 






WH(1) 


Estimated wind in target 




41-50 






WH (2) 


frame;WH(3) is usually 




51-60 






WH(3) 


zero 


2 


1-10 






BIS (1) 


Bias on range, azimuth 




11-20 






BIS (2) 


and elevation noise 




21-30 






BIS (3) 






31-40 






SIG (4) 


Initial velocity estimates 




41-50 






SIG (5) 


for filter prediction 




51-60 






SIG (6) 




3 


1-10 






SIG (1) 


Standard deviations on 




11-20 






SIG (2) 


range, azimuth, and 




21-30 






SIG (3) 


elevation noise 


4 


1-10 






SIGW(l) 


Standard deviations of 




11-20 






SIGN (2) 


assumed random forcing in 




21-30 






SI GW (3) 


x,y,and z of radar frame 


5 


1-10 






X3 (1) 


Initial position (x,y,z) 




11-20 






X3 (2) 


of aircraft in target 




21-30 






X3 ( 3) 


frame 


7 


1-10 






PHB 


Autopilot bias angle 




11-20 






TB 


Roll response time constant 




21-30 






DT 


Radar sampling interval 




31-40 






DTCON 


Control update interval 




41-50 


1 


r 


PHI 


Latitude of target 




51-60 


F10.3 


TWLD 


Time to start wild points 




61-63 


13 




MWLD 


Number of wild points 
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CARD 

8 

9 



40-49 

50 



COLUMNS 


FORMAT 


VARIABLE 


1-10 


F 1 ( 


3.3 


T3 


11-20 






THNM 


21-30 






THDM 


31-40 






THD 


41-50 






ATT 


51-60 






ATH 


1-10 






ANH 


11-20 






TUP 


21-30 






HTOL 


31-40 




t 


ATOL 


41-50 


F 1 ( 


) . 3 


S 


51-53 




13 


IB1 


54-56 




13 


IACC 


1-48 


6A8 


ATITLE- 








ZTITLE 


1-10 


Fit 


). 3 


WBLX 


11-20 






WBLY 


21-30 






BFF 


31-40 






D 




' 


t 




41-50 


F10 . 3 


W 


51-53 




13 


NTB 



REMARKS 

Nominal dive duration 

Combine to give angle 

radar and target 

Nominal dive angle 

True tang, accel. during 
dive 

Est. tang, accel. during 
dive 

Max. aircraft normal accel. 

Aircraft pitch time constant 

Angle and altitude tolerances 

for dive pullout 

Equal to 9999 

Equal to zero for level 
flight 

Equal to zero for no 
acceleration compensation 

Titles tor plots. Two 
cards per plot for two 
lines of title per plot 

Ballistic wind in x 

Ballistic wind in y 

Ballistic form factor 

Bomb diameter, inches 

Bomb weight, pounds 

Number of bombing table 
to be used 



Bombing table data as 
supplied with original 
program; this remains 
unchanged . 
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APPENDIX C: COARSE GUIDANCE INPUT DATA FORMAT DEFINITION 



CARD COLUMNS 


FORMAT 


VARIABLE 


REMARKS 


1 1 


11 




NLEG 


Number of legs 


2-10 


F9 . 


2 


RANGEI 


Range of first leg from 
radar 


11-20 


F10.4 


AZI 


Azimuth of first leg 










from radar 


21-30 






VA 


True airspeed 


31-40 






VAH 


Est. airspeed 


41-50 






VW 


True wind speed 


51-60 






VWH 


Est. wind speed 


61-70 






THETW 


True direction of wind 


71-80 


1 


r 


THETWH 


Est. direction of wind 


2 1-10 


F1C 


• 3 


PLENGTH (i) 


Length of the ith leg 


thru 11-20 

NLEG+1 






THETA(i) 


and its associated azimuth 
where i=l, NLEG 


NLEG+2 1-10 






HTG 


True initial ground heading 


11-20 






DELNM(l) 


Displacement of true position 


21-30 






DELNM(2) 


from RANGEI, AZI in x and y 


31-40 




r 


DELNM(3) 


DELNM(3) is altitude of 
aircraft in NM 


NLEG+3 1-48 


6A8 


ATITLE 


Title for plot of true. 


NLEG+4 








estimated and desired paths 


Note: All angles 


are 


in degrees from 


North, clockwise angles 


positive 


0 


- 360 degrees. All 


speeds in ft/sec. All 


lengths 


are in 


nautical miles i 


on entry. 
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